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1.  Introduction 


Dielectric  breakdown  occurs  when  a  high  voltage  or  field  is  applied  to  matter  (solid, 
liquid,  or  gas)  and  the  material  undergoes  rapid  degradation  due  to  the  stripping  of 
electrons  from  their  nuclei.  These  newly  freed  electrons  form  a  plasma  channel  as¬ 
sociated  with  high  temperatures  and  light  and  sound  emission.  In  solid  materials, 
permanent  material  failure  is  seen  along  the  breakdown  channel,  and  secondary 
fractures  may  occur  in  brittle  materials  due  to  thermal  expansion.  This  presents  a 
challenging  phenomenon  to  model  due  to  the  coupled,  multiphysics  nature  of  the 
problem.  Accurate  modeling  of  dielectric  breakdown  is  useful  in  several  applica¬ 
tions,  including  the  design  of  energy-dense  capacitors  (including  structural  capaci¬ 
tors1),  the  design  of  high-power  electronic  devices,  the  design  of  hardened  electron¬ 
ics  to  electromagnetic  pulse  attack,  and  so  forth.  Indeed,  the  multiphysics  approach 
used  here  is  especially  useful  for  structural  capacitors  as  mechanical  loading  could 
be  incorporated  into  a  simulation  along  with  electrical  loading. 

Previously,  a  dielectric  breakdown  method  was  presented  that  coupled  nonlin¬ 
ear  electro-quasi-statics,  adiabatic  heating,  and  peridynamics  using  a  finite  differ¬ 
ence  approximation  of  the  electro-quasi-static  problem  with  a  1 -phase  temperature- 
conductivity  dependency.2  Here,  the  method  is  improved  by  first  using  a  finite 
element  (FE)  discretization  of  the  electro-quasi-static  problem,  and  a  2-phase 
temperature-conductivity  material  model.  As  results  will  show,  the  simulated  break¬ 
down  patterns  better  match  experimental  breakdown  patterns  and  can  reproduce 
both  channel- like  and  tree-like  geometries. 

While  dielectric  breakdown  is  a  challenging  modeling  problem,  there  are  other  at¬ 
tempts  in  the  literature  as  well,  many  of  which  are  listed  in  Reference  2.  The  pur¬ 
pose  of  this  report  is  to  first  give  an  overview  of  the  formulation  of  the  journal 
article  cited  in  Reference  2  while  highlighting  any  differences  between  that  method 
and  the  updated  FE-based  method.  In  addition,  the  method  is  verified  against  other 
computational  methods  and  several  new  geometries  are  modeled. 

The  remainder  of  the  report  is  organized  as  follows:  First,  Section  2  details  the  for¬ 
mulation  of  the  method  including  the  physical  model,  the  discretization,  and  other 
approximations.  Next,  Section  3  verifies  the  proposed  method  against  other  compu¬ 
tational  techniques.  Section  4  presents  several  different  examples  of  the  method  on 
different  geometries,  including  3-D,  and  finally  Section  5  concludes  the  report. 
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2.  Formulation 


This  section  details  the  formulation  of  the  dielectric  breakdown  problem  in  terms 
of  3  coupled  partial  differential  equations.  The  formulation  is  similar  to  that  of  Ref¬ 
erence  2,  though  it  is  repeated  here  for  clarity.  The  main  differences  are  the  model 
used  for  the  temperature  dependence  of  the  conductivity  and  the  spatial  discretiza¬ 
tion,  which  are  discussed  in  Sections  2.3  and  2.4  respectively. 

2.1  Physical  Model 

The  model  used  here  couples  3  field  equations:  electro-quasi- statics,  the  adiabatic 
heat  equation,  and  solid  mechanics.  First,  the  electo-quasi-static  approximation 
used  has  the  form 


V>(T,  |E|)V$)  +  4v-(eV$)  =  0, 


(1) 


where  E  is  the  electric  field,  $  is  the  electrostatic  potential,  T  is  temperature,  a  is 
conductivity,  and  e  is  permittivity.  Note  that  the  equation  is  nonlinear  as  the  con¬ 
ductivity  depends  on  the  potential  (via  the  electric  field)  and  also  that  the  material 
properties  may  be  spatially  inhomogeneous  and  time-dependent.  Next,  the  temper¬ 
ature  dependence  is  defined  as 


3-T=—Q-  PHT  -  Tc), 
at  cvp 


(2) 


where  <5(-)  is  the  Dirac  delta  function,  cp  is  the  heat  capacity,  Tc  is  a  phase  change 
temperature,  (3  is  related  to  the  energy  required  for  the  phase  change,  and  Q  is  a 
heat  source  term,  which  will  be  defined  in  Section  2.2. 

Finally,  solid  mechanics  will  be  used  to  model  forces  on  a  solid  body  due  to  thermal 
expansion  and  electrostatic  forces  (Lorentz  and  Kelvin),  though  peridynamics  will 
be  used  rather  than  standard  elastodynamics  as  we  wish  to  model  fracture  and  fail¬ 
ure  of  a  material  due  to  high  strains  and  temperatures.  Peridynamics  is  a  nonlocal 
formulation  of  elastodynamics  that  naturally  incorporates  discontinuities  that  arise 
from  fracture.3,4  Bond-based  peridynamics  may  be  stated  as 


(3) 
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where  the  so-called  microforce  function  is  given  as 


f(77,  tT)=c  [s(rj,  i)  -  a  AT]  (4) 

l*7  +  CI 

where  AT  is  the  temperature  difference  relative  to  ambient  Tam b  and  a  is  the 
isotropic  thermal  expansion  coefficient.  The  peridynamic  stretch  s  is  defined  as 


s(»7,£) 


1*7  +  11 

— 

ICI 

1C 

(5) 


where 

^  =  x'  -  x, 


and 


T]  =  ll'  —  U, 


(6) 

(7) 


and  constant  c  is  defined  as 


6  E 

7 r<53  (1  —  u)  ’ 


(8) 


for  2-D  plane  stress,  or 

677 

c  =  7 r54  (1  -  2u)  ’  (9) 

for  3-D  in  terms  of  Young’s  modulus  E  and  Poisson’s  ratio  v.  The  only  difference 
between  this  formulation  and  that  in  Reference  4  is  that  isotropic  thermal  expansion 
is  included. 


Damage  is  modeled  in  bond-based  peridynamics  with  a  bond-breaking  scheme, 
wherein  the  microforce  is  set  to  zero  if  a  pair  of  points  has  ever  had  a  stretch  that 
exceeds  a  given  critical  value.  The  microforce  is  then  modified  as 


f(77,  £,  T)  =  ch{t ,  0  [s(rj,  0  -  aAT]  (10) 

where  the  health  h  of  a  given  bond  may  be  expressed  as 


if  ~ 

otherwise 


aAT(f')  <  s0  and  Tavg(f')  <  Tc  for  all  0  <  t'  <  t 

(11) 
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and  the  critical  stretch  s0  is 


for  2-D  plane  stress  or 


so  — 


(12) 

(13) 


for  3-D.  In  Eq.  11,  excessive  temperature  will  also  fail  a  bond,  and  Tavg  is  the  av¬ 
erage  temperature  of  the  2  nodes  x  and  x'  and  Tc  is  a  given  critical  temperature. 
Finally,  peridynamics  has  no  concept  of  discrete  cracks  or  fracture  surfaces,  though 
damage  can  be  defined  as  the  ratio 


d 


In,  h  dV> 


(14) 


This  value  is  useful  for  postprocessing  and  will  also  be  used  in  a  modified  permit¬ 
tivity  as  described  in  the  next  section. 


2.2  Coupling 

The  main  field  equations  (Eqs.  1,  2,  and  3)  are  coupled  together  in  several  ways. 
First,  electrostatic  forces  couple  the  electric  field  to  the  displacement,  the  first  of 
which  is  the  Lorentz  force,  defined  as 


FL  =  gE  =  V-(eE)E,  (15) 

where  q  is  the  charge  density.  This  force  acts  on  free  charges  in  the  material  in 
areas  of  high  conductivity.  The  second  force  is  the  Kelvin  force,  or  the  force  on  a 
dielectric,  defined  as 


Fk  =  P-  VE=  (e-eo)E-  VE,  (16) 

where  P  is  the  polarization  and  e0  is  the  permittivity  of  free  space.  The  equation  of 
motion  is  then  modified  to 

q2 

—  Fpd  +  FL  +  Fk,  (17) 

at 2 

where  FPD  is  the  peridynamic  force  defined  on  the  right-handside  of  Eq.  3. 
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Here,  the  peridynamic  damage  is  coupled  to  the  permittivity  by  the  linear 
relationship 

e  =  eo  [er  (1  —  d)  +  d] .  (18) 

The  intent  is  to  model  the  effect  of  damage  in  a  material  as  the  formation  of  voids 
so  that  a  fully  damaged  material  point  (with  d  —  1)  acts  as  free  space. 

The  heat  source  term  Q  in  Eq.  2  is  derived  from  Joule  heating,  or  the  heat  generated 
from  current  flowing  through  a  dissipative  material,  and  is  given  by 

Q  =  J-E  =  a|E|2.  (19) 


Finally,  thermal  expansion  (defined  in  Eq.  4)  couples  the  temperature  to  the  me¬ 
chanical  displacement. 

2.3  Conductivity  Model 

As  mentioned  above,  the  conductivity  is  dependent  on  both  the  electric  field  E  and 
the  temperature  T.  The  model  used  here  is  exponential  and  given  by 

a(T,\E\)  =  a0f(T)e^,  (20) 


where  a0  is  a  base  conductivity  with  zero  applied  field,  and  /  is  a  function  de¬ 
scribing  the  temperature  dependence.5  The  temperature  dependence  is  a  2-phase 
Arrhenius  model  given  by6 


f(T ) 


Aie~Bl/T  T  <TC 
A2e~B 2/t  T>Tc  ’ 


(21) 


where  the  A,;  and  B,  are  constants  of  the  model.  The  constants  A\  and  B\  are  as¬ 
sociated  with  the  weakly  conductive  phase  and  A2  and  B2  are  associated  with  the 
highly  conductive  phase,  meaning  that  A±  A2. 


2.4  Discretization 


The  spatial  discretization  of  Eq.  3  is  a  basic  collocation  method  as  is  used  perva- 
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sively  in  the  literature,4  given  by 


N i 


Fp„=  V 


j\rje'Hri 


(22) 


where  V3  is  the  volume  of  node  j.  The  temporal  discretization  is  a  Velocity  Verlet 
method,  given  as 


\k+1/  2  =  \k  +  — afc, 

u k+1  =  uk  +  Afvfc+1/2, 
„fc+i  _  Fpd  +  Fk  +  Fl 

d  —  j 

P 

yfc+l  _  yfe+i/  2  _|_ 


(23) 


The  main  difference  between  this  work  and  that  of  Reference  2  is  that  the  electro- 
quasi-static  equation  is  solved  using  a  finite  element  method  (FEM).  First,  the  FE 
formulation  is  given  as 

Y  4>n  f  ° V tm  •  VbndV  -  Y  °"  j  '  'VKdV  =  0,  (24) 

where  the  bn  and  trn  are  basis  and  testing  functions  respectively,  Qm  is  the  support 
of  the  mth  testing  function,  and  <i>n  are  the  nodal  potentials.  As  usual,  the  region 
is  broken  into  elements,  where  here  quadrilateral  (2-D)  and  hexahedral  (3-D)  ele¬ 
ments  are  used  with  bilinear  and  trilinear  basis  and  testing  functions.  Further,  the 
peridynamic  nodes  are  located  at  the  centroids  of  the  elements. 

An  entry  in  a  system  matrix  is  defined  as 

[D a\mn  =  j  aVtm  ■  VbndV,  (25) 

with  a  being  a  dummy  variable  representing  either  a  or  e.  The  system  matrices  are 
formed  in  the  usual  way  using  assembly  of  element  matrices. 

Next,  the  time  dependence  of  the  adiabatic  heating  model  is  discretized  using  ex- 
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plicit  forward  Euler  as 


Tk+ 1  =  Tk  +  — Qk  -  A t/35  (Tk  -  Tc)  .  (26) 

cPP 

To  facilitate  the  evaluation  of  Eq.  26,  the  delta  function  is  regularized  with  a  Gaus¬ 
sian  as  follows 

m  «  (27) 

yjna 

where  a  is  the  variance  and  is  a  parameter  of  the  model. 

2.5  Linearization 

The  electro-quasi-static  equation  is  nonlinear  in  the  potential  because  the  conduc¬ 
tivity  depends  on  the  magnitude  of  the  electric  field.  This  equation  may  be  solved  in 
a  few  ways,  and  here  we  use  a  linearization  procedure  along  with  a  backward  Euler 
temporal  discretization.  First,  the  FE-discretized  approximation  may  be  expressed 
as 

=  0,  (28) 

where  system  matrices  Da  and  Df  are  defined  in  Eq.  25  and  $  is  a  vector  of  nodal 
unknowns  representing  the  potential.  Now,  the  temporal  derivative  is  applied  to  both 
the  potential  and  the  permittivity,  as  the  permittivity  is  time-dependent  (because  it 
depends  on  peridynamic  damage),  giving 

+  Dc,$  +  De4$  =  0,  (29) 

at 

where  e'  is  the  temporal  derivative  of  the  permittivity.  Backward  Euler  is  used  to 
discretize  the  temporal  dependence,  giving 


(AtD^  +  2Defc  -  Defc-i)  =  Defc<E»fc_1.  (30) 

Note  that  superscript  k  indicates  the  variable  at  the  kth  time  step  and  that  the  dis¬ 
cretization  is  also  applied  to  the  time  derivative  on  the  permittivity. 

Now,  Eq.  30  may  be  linearized  by  expressing  the  conductivity  times  the  electric 
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field  as  a  Taylor  series  about  the  field  at  the  previous  time  step,  which  gives 


a(Tk,  |Efc|)Efc  «  a0f(Tk ) 


k\ei\vk 


Ek 


E' 


k~l  gfc-1 


7- 


E 


3- - (Efc  -  Efc_1) 


Finally,  this  linearization  is  inserted  into  Eq.  30  giving, 


A#D  k  +  AtD-k  +  2D fk  —  D./c-i  )  = 


Defc  +  AtDak  )  <f>k~\ 


(31) 


(32) 


where  D„/.  is  the  FE  matrix  associated  with  the  first  term  on  the  right-handside  of 
Eq.  31  and  D is  the  FE  matrix  associated  with  the  second  term.  Note  that  the 
second  term  is  a  tensor,  and  so  the  definition  of  Eq.  25  is  extended  to 


[D, 


Vim A  [ Vbn}TdV . 


(33) 


2.6  Nonlocal  Force  Computation 

The  final  step  in  the  formulation  is  the  computation  of  the  electrostatic  forces  FK 
and  Fl.  Unfortunately,  they  require  that  the  basis  functions  be  at  least  twice  differ¬ 
entiable  (see  Eq.  15  and  Eq.  16)  though  here  only  bilinear  (or  trilinear  in  3-D)  basis 
functions  are  used.  Rather  than  using  a  higher-order  basis,  a  nonlocal  approxima¬ 
tion  to  the  force  computation  will  be  used.  The  following  is  akin  to  the  state -based 
peridynamic  formulations7  and  concepts  from  nonlocal  calculus.8 

The  first  step  in  computing  the  electrostatic  forces  is  to  compute  the  electric  field  E, 
which  may  be  obtained  by  using  the  derivative  of  the  basis  functions.  Here,  the  field 
will  be  evaluated  at  the  centroid  of  each  element,  coincident  with  the  peridynamic 
nodes.  Next,  the  gradient  of  the  field  is  needed  for  computing  the  Kelvin  force  and 
is  approximated  as 


VE 


/  c  (|£|)  (E7  —  E)  ®  £dx' 
J-H* 


K 


-i 


(34) 


where  c  (|£|)  is  a  shape  function  (which  is  simply  constant  over  the  horizon  here) 
and  the  same  horizon  is  used  as  the  peridynamic  method  discussed  above,  and  K  is 
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a  shape  tensor  given  by 


c(|£|)£<8>  £dx'. 


(35) 


Finally,  the  divergence  of  the  electric  flux  density,  D 
the  Lorentz  force,  and  is  given  as 


eE,  is  needed  to  compute 


V  •  D  trace 


[  c(l£l)(D' 


D)  <g)  £d,x' 


(36) 


2.7  Algorithm 

To  summarize,  the  previously  discussed  methods  are  assembled  in  the  following 
way: 

1 .  Specify  any  initial  electrical  material  properties  given  by  the  problem  geom¬ 
etry  and  compute  initial  potential 

2.  Update  the  displacement  and  velocity  based  on  acceleration  using  Eq.  23 

3.  Update  the  temperature  using  Eq.  26 

4.  Update  the  bond  health  due  to  stretch  and  temperature  (Eq.  11) 

5.  Update  the  permittivity  based  on  the  damage  using  Eq.  18 

6.  Update  the  conductivity  based  on  the  electric  field  and  temperature  using 
Eq.  31 

7.  Compute  the  electrostatic  potential  (<l>)  from  Eq.  32,  with  specified  boundary 
conditions 

8.  Compute  the  electrostatic  forces  from  the  potential  using  nonlocal  approxi¬ 
mations  (Fl  and  Fk)  with  Eq.  15  and  Eq.  16 

9.  Update  the  velocity  and  acceleration  using  all  forces  FPD,  FL,  and  FK 

10.  Repeat  starting  at  Item  2 

This  method  uses  a  number  of  approximations  in  not  only  the  spatial  and  temporal 
discretizations,  but  also  in  the  computation  of  the  electrostatic  forces  and  other 
regularizations.  It  is  therefore  important  to  verify  various  aspects  of  the  method,  as 
discussed  in  the  next  section. 
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3.  Verification 


A  few  approximations  used  in  the  approach  were  verified  using  other  computa¬ 
tional  methods.  First,  the  computation  of  the  Kelvin  force  was  verified  against  a 
Mathematica  solution  wherein  a  fictitious  force  was  specified  and  the  permittivity 
and  potential  that  generate  that  force  were  computed.  Next,  the  linearization  of  the 
nonlinear  conductivity  model  was  verified  with  a  fixed-point  iteration  scheme. 

3.1  Nonlocal  Kelvin  Force  Computation 

The  computation  of  the  Kelvin  force  uses  a  nonlocal  formulation  due  to  the  use  of 
bilinear  basis  functions,  which  is  verified  here  using  the  following  approach:  First,  a 
fictitious  force  is  specified  in  1-D.  Next,  the  permittivity  and  potential  that  generate 
that  force  are  solved  numerically.  Finally,  the  permittivity  and  potential  are  extended 
to  2-D  (constant  along  one  dimension)  and  the  Kelvin  force  is  computed  using 
the  nonlocal  approach.  The  resulting  force  can  then  be  compared  to  the  original 
fictitious  force. 

The  force  used  here  has  the  following  form 

fK(x)  =  -10_11sin7nr,  (37) 

with  0  <  x  <  1  and  the  permittivity  e  and  potential  0  are  defined  via  the  coupled 
ordinary  differential  equations  (ODEs) 

e  (x)  0"  ( x )  +  e'  (x)  (f>  (x)  =  0  (38) 

and 

e0  [e  (x)  -  1]  0'  (x)  0"  (x)  +  /K  (a?)  =  0,  (39) 

subject  to  the  boundary  conditions 


<M0)  =  0, 

0(1)  =  1,  (40) 

e(0)  =  2. 

Equations  38  and  39  were  solved  for  permittivity  e  and  potential  0  numerically  us¬ 
ing  Mathematica’s  built-in  ODE  solver  and  are  shown  in  Fig.  1.  As  stated  above, 
these  solutions  were  extended  to  2-D  along  one  dimension  and  used  with  the  non- 
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local  Kelvin  force  computation. 


(a)  (b) 

Fig.  1  Permittivity  a)  and  potential  b)  that  generate  the  Kelvin  force  of  Eq.  37 


Figure  2  shows  the  nonlocal  approximation  of  the  Kelvin  force  in  2-D  using  100 
elements  along  each  axis  and  a  horizon  size  of  3 Ax.  To  compare  against  the  1-D 
solution,  a  slice  along  the  center  is  used  as  shown  in  Fig.  3.  Figure  3  also  shows  the 
Kelvin  force  as  computed  by  a  finite  difference  implementation.2  The  error  conver¬ 
gence  of  the  Kelvin  force  computation  was  also  assessed  as  shown  in  Fig.  4,  where 
both  the  finite  difference  and  FE  methods  are  shown.  As  can  be  seen,  the  error  in 
the  Kelvin  force  decreases  with  an  increasing  number  of  elements. 

3.2  Linearization 

The  linearization  described  above  was  verified  against  a  fixed-point  iteration  solu¬ 
tion  of  the  system  of  equations  given  in  Eq.  30.  To  perform  the  fixed-point  iteration, 
Eq.  30  is  first  rewritten  as 

F($)$=b,  (41) 

where  F  is  the  matrix  defined  on  the  left-handside  of  Eq.  30  and  is  a  function  of  <F. 
and  b  is  the  right-handside  of  Eq.  30.  This  leads  to  the  fixed-point  iteration 

$,+1  =  F($I)-1b,  (42) 

and  superscript  l  refers  to  a  step  in  the  fixed-point  iteration,  not  the  time  step.  This 
iteration  was  found  to  converge  if  the  problem  was  well-behaved  (i.e.,  at  low  enough 
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voltages  so  that  no  damage  occurred). 


Fig.  2  Nonlocal  approximation  of  Kelvin  force  in  2-D 
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Fig.  3  Comparison  of  Kelvin  force  computed  using  FEM  and  finite  differences 


For  comparison,  the  fixed-point  iteration  solution  was  used  rather  than  the  lineariza¬ 
tion  scheme  at  a  timestep  size  of  125  ps  for  the  point-plane  geometry  shown  in 
Fig.  5.  Only  the  blue  region  (dielectric)  of  Fig.  5  was  meshed  and  an  average  ele¬ 
ment  edge  length  of  50  pm  was  used.  The  gold  regions  represent  perfect  conduc¬ 
tors  and  are  boundary  conditions  in  the  model.  A  homogeneous  Neumann  boundary 
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condition  was  used  on  the  dielectric-free  space  interfaces  (top,  left,  and  right  edges). 
In  addition,  the  applied  voltage  had  the  form 

$applied(£)  =  V"raax  (l  -  e~t/T )  ,  (43) 

where  time  constant  r  is  0.3  ps  unless  specified  otherwise. 


Fig.  4  Error  convergence  of  Kelvin  force  computed  using  FEM  and  finite  differences 


The  linearized  version  was  then  run  at  3  discretizations  to  a  final  time  of  100  ns:  1 
ns,  500  ps,  and  250  ps.  At  the  time  step  corresponding  to  100  ns  (800  for  the  fixed- 
point  solution  and  100,  200,  or  400  for  the  linearization)  the  error  in  the  electric 
field  distribution  was  measured  relative  to  the  fixed-point  iteration  solution.  Figure  6 
shows  the  convergence  of  the  linearization  to  the  fixed-point  solution  for  the  3  time- 
step  sizes.  As  can  be  seen  from  the  figure,  decreasing  the  time-step  size  improves 
the  accuracy  of  the  linearized  solution.  Finally,  Fig.  7  shows  the  log  error  in  E  of 
the  linearized  solution  with  a  1  ns  time-step  size  after  100  ns.  As  expected,  the  error 
is  highest  near  the  highest  field  concentration,  which  is  where  the  conductivity  will 
have  the  highest  variation. 
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cK7.5mm 


m 


/7=10mm 


w=10mm 


Fig.  5  Schematic  of  the  point-plane  geometry 


#  Time  steps 


Fig.  6  Convergence  of  the  linearization  of  Eq.  30  to  a  fixed-point  iteration  solution 


4.  Results 

As  a  first  set  of  tests,  point-plane  geometries — wherein  a  point  probe  is  embedded 
in  a  dielectric  and  placed  above  a  ground  plane — were  modeled.  First,  a  flat-tipped 
electrode  was  modeled  twice,  first  with  material  properties  that  generate  a  straight 
channel-like  breakdown  pattern,  and  second  with  properties  that  generate  a  tree- 
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like  structure.  Next,  a  sharply  tipped  electrode  is  modeled.  The  last  2-D  example  is 
a  parallel  plate  capacitor  with  a  composite  dielectric  and  a  small  conductive  flaw. 
Finally,  a  3-D  point-plane  geometry  is  tested  with  a  cylindrical  dielectric. 


x  10 3 


Fig.  7  Error  in  E  of  the  linearization  with  a  time-step  size  of  1  ns  vs.  a  fixed-point  solution 


Unless  specified  otherwise,  all  material  properties  used  in  this  section  are  given  in 
Table  1  and  the  applied  voltage  has  the  form  given  in  Eq.  43. 

4.1  Flat  Electrode 

This  example  repeats  the  results  from  Reference  2,  though  now  with  the  2-phase 
model  described  previously  (a  schematic  of  the  geometry  is  shown  in  Fig.  5).  Aside 
from  the  2-phase  temperature-conductivity  relation,  the  material  properties  and  ge¬ 
ometry  are  identical  to  those  described  in  Reference  2.  A  mesh  of  quadrilateral 
elements  is  used,  with  average  edge  length  of  50  pm. 

The  results  are  shown  in  the  following  set  of  figures:  Fig.  8  shows  the  peridynamic 
damage  at  4  time  steps  between  4  and  5.5  ps  and  Fig.  9  shows  the  temperature  and 
conductivity  at  5  ps.  Note  that  the  conductivity  is  shown  on  a  logarithmic  scale. 

In  contrast  to  the  results  in  Reference  2,  a  more  clearly  defined  breakdown  chan¬ 
nel  is  formed,  most  likely  due  to  the  2-phase  model  used  for  the  temperature- 
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conductivity  dependence.  Perpendicular  fractures  are  also  seen,  which  are  due  to 
the  thermal  expansion  in  the  breakdown  channel  applying  a  force  to  the  surround¬ 
ing  material.  This  effect  is  seen  in  dielectric  breakdown  experiments.9 


Damage 

:1.000e+00 


E0.75 


E0.5 


z0.25 

ur0.000e+00 


(d) 


Fig.  8  Peridynamic  damage  in  the  flat-tipped  point-plane  model  at  4  ps  a),  4.5  ps  b),  5  ps  c), 
and  5.5  ps  d) 


4.2  Tree-Like  Breakdown 

The  same  geometry  was  used  as  in  the  previous  example,  though  now  the  follow¬ 
ing  material  properties  were  altered:  The  electric  field-conductivity  coupling  coef¬ 
ficient  was  set  to  7  =  0,  the  ambient  temperature  was  increased  to  Tamb  =  670  K, 
the  maximum  voltage  was  Vo  =  50  kV,  and  the  initial  conductivity  was  raised  to 
cr0  =  1  S/m.  Note  that  now  the  electro-quasi-static  equation  is  linear.  These  mate¬ 
rial  properties  lead  to  the  tree-like  breakdown  pattern  seen  in  Fig.  10,  which  shows 
the  peridynamic  damage  at  1.76  ps.  It  is  hypothesized  that  the  tree-like  pattern  (ver¬ 
sus  the  straighter,  channel-like  pattern  seen  in  Fig.  8)  is  because  the  nonlinearity  in 
the  previous  example  generates  a  channel  of  high  conductivity,  and  breakdown  can 
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only  occur  in  this  channel.  In  this  example,  breakdown  may  occur  anywhere  due  to 
the  high  conductivity  throughout  the  material  and  so  a  tree-like  pattern  develops. 
Further,  the  stochastic,  unstable  damage  pattern  is  most  likely  due  to  the  irregular 
mesh. 


(a)  (b) 

Fig.  9  Temperature  a)  and  conductivity  b)  in  the  flat-tipped  point-plane  model  at  5  ps 


4.3  Sharply  Pointed  Electrode 

Here,  nearly  the  same  model  as  that  discussed  in  Section  4.2  is  presented,  though 
now  the  tip  of  the  electrode  is  sharp  rather  than  flat.  Overall,  the  results  are  nearly 
the  same,  though  the  breakdown  channel  is  thinner  at  first.  In  addition,  breakdown 
occurs  over  a  longer  period  of  time  versus  the  flat-tipped  model.  The  peridynamic 
damage  is  shown  in  Fig.  11  between  2.5  and  10  ps. 

The  same  input  parameters  from  the  tree-like  breakdown  simulation  discussed  in  the 
previous  section  were  repeated,  resulting  in  the  damage  pattern  shown  in  Fig.  12. 
Again,  a  tree-like  pattern  is  seen,  though  now  only  one  branch  originates  from  the 
electrode  tip.  As  before,  it  is  believed  that  the  unstable,  tree-like  damage  pattern  is 
due  to  the  irregular  mesh. 
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Table  1  Physical  constants 


Constant 

Value 

Units 

Mechanical 

P 

2400 

kg  m  3 

E 

72 

GPa 

8 

0.15 

mm 

G0 

5 

J  m-2 

Electrical 

eo 

8.85  x  10"12 

F  m_1 

er 

20 

- 

CTO 

io-19 

Sm_1 

Thermal 

Cp 

800 

Jkg^K"1 

Tc 

300 

K 

Tirnb 

1000 

IC 

Coupled 

a 

9  x  10~6 

K"1 

7 

5  x  10"8 

mV-1 

a 

5 

K 

P 

8  x  107 

K2 

A! 

30 

- 

^2 

3  x  104 

- 

B\,  B2 

1.2  x  103 

K 

4.4  Composite  Capacitor  with  Conductive  Flaw 

The  final  2-D  example  is  a  parallel-plate  capacitor  with  a  composite  dielectric.  The 
model  used  2  dielectric  materials  along  with  a  small,  conductive  flaw  in  the  inner 
dielectric,  as  shown  in  Fig.  13,  and  the  conductors  were  located  along  the  top  and 
bottom  edges.  The  model  is  15  mm  by  10  mm  with  a  3-mm-radius  circular  dielectric 
located  at  the  center.  The  outer  region  had  a  relative  permittivity  of  10  and  a  thermal 
expansion  coefficient  of  9  x  10-6  and  the  inner  region  had  a  relative  permittivity  of 
20  and  a  thermal  expansion  coefficient  of  20  x  10-6.  Within  the  inner  dielectric  is 
a  circular,  conductive  inclusion  (conductivity  10”3  S/m)  with  a  radius  of  75  pm.  In 
addition  to  having  different  dielectric  constants,  the  2  materials  had  different  ther¬ 
mal  expansion  coefficients  to  induce  fracture  at  the  interface.  The  maximum  voltage 
used  in  this  simulation  was  4  MV  with  the  same  time  constant  used  previously. 
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Fig.  10  Tree-like  breakdown  pattern  for  linear  breakdown  problem 


In  this  model,  breakdown  initiates  at  the  inclusion,  since  the  inclusion  concentrates 
the  electric  field  (as  shown  in  Fig.  14)  at  a  small  point.  The  breakdown  then  prop¬ 
agates  towards  the  boundaries  as  shown  in  Fig.  15.  This  damage  is  associated  with 
high  temperatures  due  to  breakdown,  though  there  are  other  fractures  induced  by 
thermal  expansion  visible  especially  around  the  inner  dielectric.  Finally,  the  temper¬ 
ature  and  conductivity  are  shown  in  Fig.  16  and  Fig.  17.  Note  that  the  conductivity 
is  shown  on  a  log  scale. 

4.5  Point-Plane  Geometry  in  3-D 

A  3-D  model  was  generated  (shown  in  Fig.  1 8),  which  consisted  of  a  sharply  pointed 
probe  of  length  1  mm  and  radius  0. 1  mm  embedded  in  a  dielectric  cylinder  with  a 
height  of  4  mm  and  a  radius  of  2  mm.  As  before,  the  voltage  is  applied  to  the  probe, 
while  the  entire  bottom  surface  of  the  cylinder  is  the  ground  plane.  While  the  ma¬ 
terial  properties  were  the  same  as  above,  the  temporal  discretization  was  reduced 
to  At  =  0.75  ns  and  the  maximum  voltage  was  reduced  to  1.1  MV.  The  results 
show  a  thin  breakdown  channel  forming  at  approximately  375  ns  and  propagating 
more  slowly  than  the  2-D  results.  First,  the  electric-field  magnitude  and  conductiv¬ 
ity  (on  a  logarithmic  scale)  are  shown  in  a  cutaway  view  of  the  dielectric  at  2.2  ps 
in  Fig.  19.  Next,  Figure  20  shows  a  cutaway  view  of  the  peridynamic  damage  in 
the  center  of  the  cylinder  at  4  instances  of  time  between  2.2  ps  and  4.5  ps.  Compar- 
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ing  Fig.  19b  with  Fig.  20a  shows  the  peridynamic  damage  that  is  associated  with 
breakdown,  and  not  other  effects  such  as  thermal  expansion.  Further,  a  threshold 
view  of  the  peridynamic  damage  is  shown  in  Fig.  21  at  4.5  ps,  first  in  a)  with  a 
threshold  on  the  damage  of  0.25  and  second  in  b)  with  a  threshold  on  temperature 
(1000  K).  In  other  words,  Fig.  21b  shows  the  damage  in  regions  with  high  tempera¬ 
ture,  above  the  specified  breakdown  temperature  threshold.  The  difference  between 
the  2  figures  illustrates  the  different  types  of  damage:  The  thin  channel  in  Fig.  21a 
is  due  to  high  temperature  (resulting  from  high  conductivity  as  shown  in  Fig.  21b) 
and  the  fractures  at  the  bottom  of  the  cylinder  are  due  to  thermal  expansion  as  the 
temperature  in  that  region  is  below  the  breakdown  threshold. 


(c)  (d) 

Fig.  11  Peridynamic  damage  in  the  sharp-tipped  point-plane  model  at  2.5  ps  a),  5  ps  b),  7.5  ps 
c),  and  10  ps  d) 


Finally,  a  tree-like  breakdown  pattern  was  generated  in  3-D  by  first  using  a  flat- 
tipped  electrode  embedded  in  a  dielectric  cylinder  with  the  dimensions  the  same 
as  those  stated  previously.  In  this  case,  the  spatial  discretization  size  was  reduced 
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to  an  average  element-edge  length  of  40  |am.  The  material  properties  and  initial 
temperature  were  the  same  as  those  in  the  tree-like  breakdown  example  in  2-D, 
though  here  a  maximum  voltage  of  10  kV  was  used.  Figure  22  shows  the  damage 
at  3.3  ps  by  superimposing  a  threshold  view  of  the  damage  (with  a  threshold  of 
0.35)  on  a  cut-away  view  of  the  damage.  As  before,  a  tree-like  breakdown  pattern 
is  evident  that  initiates  at  the  flat-tipped  electrode. 


Damage 


1 .000e+00 


0.75 

0.5 

0.25 


0.000e+00 


Fig.  12  Tree-like  breakdown  pattern  for  pointed  electrode  at  1.5  ps 


5.  Conclusions 

A  multiphysics,  hybrid  FE-peridynamics  method  was  presented  for  solving  dielec¬ 
tric  breakdown  problems.  The  method  coupled  3  field  equations:  electro-quasi- 
statics  for  solution  of  the  electrostatic  potential,  adiabatic  heating  for  solution  of 
the  temperature,  and  peridynamics  for  displacement.  These  3  equations  were  cou¬ 
pled  in  various  ways:  Lorentz  and  Kelvin  electro-static  forces,  a  temperature-  and 
electric-field-dependent  conductivity  model,  damage-dependent  permittivity,  and 
Joule  heating.  Results  were  presented  for  different  geometries  in  2-D  and  3-D  and 
some  characteristics  were  consistent  with  experiments.  The  most  encouraging  result 
is  that  the  method  is  capable  of  generating  both  channel-like  and  tree-like  break¬ 
down  geometries  depending  on  the  material  properties  chosen.  The  channel-like 
solutions  were  generated  using  the  nonlinear  electric-field-conductivity  dependence 
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and  low  initial  conductivity  and  temperature.  For  the  point-probe  geometries,  it  ap¬ 
pears  the  field  configuration  raises  the  conductivity  along  a  straight  path  between 
the  electrode  and  ground  plane,  thus  only  allowing  breakdown  in  this  area.  The 
tree-like  solutions  were  generated  with  a  constant  conductivity  and  high  initial  con¬ 
ductivity  and  temperature.  In  this  case,  the  higher  background  conductivity  allows 
for  breakdown  to  occur  on  different,  unstable  paths.  These  results  compare  well  to 
Lichtenberg  figures,  some  of  which  are  generated  by  bombarding  a  polymer  plate 
with  electrons  and  subsequently  applying  a  large  voltage.  The  embedded  electrons 
raise  the  charge  in  the  material,  which  is  similar  to  raising  its  conductivity. 


Fig.  13  Composite  parallel  plate  capacitor  model  with  mesh 


While  this  approach  has  generated  results  that  resemble  experiments,  several  im¬ 
provements  could  be  made.  First,  the  displacement  is  not  coupled  to  the  electro¬ 
static  problem,  meaning  the  deformation  of  the  body  does  not  affect  the  solution  of 
the  electro- static  problem.  In  the  materials  used  here,  the  displacements  are  small 
enough  to  be  neglected,  though  for  different  materials,  such  as  electro-active  poly¬ 
mers,  large  deformations  would  necessitate  more  accurate  handling  of  this  cou¬ 
pling.  Second,  only  a  2-phase  model  is  used  for  the  temperature-conductivity  de¬ 
pendence,  representing  low  and  high  conductivity.  In  reality,  a  solid  undergoing 
dielectric  breakdown  may  transition  through  all  phases  of  matter,  from  solid  to  liq- 
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uid,  gas,  and  plasma.  Future  implementations  may  incorporate  material  models  for 
each  phase,  not  just  in  the  temperature-conductivity  dependence,  but  also  in  the 
equation  of  motion  and  elastic  and  thermal  material  properties.  Third,  the  nonlinear 
electric-field-conductivity  dependence  is  based  on  experimental  support  that  may 
not  cover  the  entire  range  of  electric-field  magnitudes  used  here,  and  so  in  the  fu¬ 
ture  more  sophisticated  field-dependent  conductivity  models  should  be  explored. 
Finally,  magneto- static  forces  may  also  be  included,  though  it  remains  to  be  seen 
whether  or  not  the  magnitude  of  such  forces  would  impact  the  solution. 


Fig.  14  Magnitude  of  the  electric  field  for  the  composite  capacitor  model  at  0.5  ps 
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Fig.  15  Peridynamic  damage  in  the  composite  capacitor  model  at  5  ps  a),  5.5  ps  b),  and  6  ps  c) 
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(a) 


(b) 


Fig.  16  Temperature  in  the  composite  capacitor  at  5.5  ps  a)  and  5.5  ps  b) 
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(a) 


Conductivity 


(b) 


Fig.  17  Conductivity  in  the  composite  capacitor  at  5.5  ps  a)  and  5.5  ps  b) 
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Fig.  18  Cylindrical  point-plane  model  (cutaway) 


(a) 


(b) 


Fig.  19  Electric-field  magnitude  a)  and  conductivity  b)  in  the  3-D  point-plane  model  at  2.2  ps 
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(d) 


Fig.  20  Peridynamic  damage  in  the  3-D  point-plane  model  at  2.25  ps  a),  3  ps  b),  3.75  ps  c),  and 
4.5  ps  d) 
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(a)  (b) 


Fig.  21  Peridynamic  damage  in  the  3-D  point-plane  model  at  4.5  ps  shown  with  a  threshold 
on  damage  a)  and  temperature  b) 
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Fig.  22  Damage  for  the  flat-tipped  probe  in  3-D  at  3.3  ps 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


1- D  1 -dimensional 

2- D  2-dimensional 

3- D  3-dimensional 
FE  finite  element 

FEM  finite  element  method 
ODE  ordinary  differential  equation 
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